Initialisation

## [1] "number of cores available = 1"
#Phi[1] ; eta = valeur de fin
#Phi[2] = valeur du noeud
#Phi[3] = echelle
m <- function(t, eta, phi) (phi[,1] + eta)/(1+exp((phi[,2]-t)/phi[,3]))
#=======================================#
param <- list(sigma2 = 0.05,
              rho2 = 0.1,
              mu = c(5,4,0.5),
              omega2 = c(0.5,0.1,0.01))

F. <- length(param$omega2) #dimension de phi
#=======================================#
t <- seq(2,6, length.out = 10) #value of times

# --- longitudinal data --- #
G = 40 ; ng = 8 #nombre de groupe et d'individu par groupe
n <- G*ng*length(t) ; N <- G*ng  #nombre total de data, et #nombre total d'individu

dt_NLME <- NLME_obs(G, ng, t, param, m)
Y <- dt_NLME$obs

SAEM avec simulation par MCMC

\[\log f(Y, Z ; \theta) = \langle \Phi(\theta) ; S(Y,Z) \rangle - \psi(\theta)\]

#Petite fonction pour retourner rapidement l'appel Phi[attr(Phi, i)] où i est '1', '2', ...,
`%a%` <- function(x,var){
  if(length(var)== 1) return(x[ attr(x,as.character(var)) ])
  lapply(var, function(v) x%a%v) %>% unlist
}


Phi <- fct_vector(function(sigma2, rho2, mu, omega2) -n/(2*sigma2),
                function(sigma2, rho2, mu, omega2) -N/(2*rho2),
                function(sigma2, rho2, mu, omega2) -G/(2*omega2),
                function(sigma2, rho2, mu, omega2) G*mu/omega2, dim = c(1,1,F., F.) )

S <- fct_vector(function(eta, phi) mean((Y - get_obs(m, dt_NLME, eta = eta, phi = phi) )^2 ),
              function(eta, phi) mean(eta^2),
              function(eta, phi) apply(phi^2, 2, mean),
              function(eta, phi) apply(phi  , 2, mean), dim = c(1,1,F., F.) )

Metropolis Hastings

loglik.phi <- function(x, eta, Phi)
{
  id <- c(1,3,4) #indice de S_1 et S_{3,.} puis S_{4,.}
  sum( Phi%a%id * S(eta, x, i = id) )
}
loglik.eta <- function(x, phi, Phi) 
{ 
  id <- c(1,2)
  sum( Phi%a%id * S(x, phi, i = id) )
}

SAEM

Initialisation

M <- 1 #nombre de simulation

# ---  Initialisation des paramètres --- #
para <- param %>% lapply(function(x) x + runif(1))
#para$rho2 = 0.2 ; para$mu <- c(6,3,1) ; para$omega2 <- rep(.1,3)

# --- Initialisation des chaines MC : Z_0 ---
Z <- list(eta = 1:M %>% lapply(function(i) rnorm(G*ng, 0, para$rho2)  %>% matrix(ncol = 1) ),
          phi = 1:M %>% lapply(function(i) matrix(rnorm(F.*G, para$mu, para$omega2), nrow = F.) %>% t ) )

Etape simulation et maximisation du SEAM

sim <- function(niter, Phih, eta, phi)
{
  M <- length(phi)

  eta <- 1:M %>% lapply( function(i)
    MH_High_Dim_para_future(niter, eta[[i]], sd = .036,                 loglik.eta, phi[[i]], Phih, cores = ncores-1 ))
  
  phi <- 1:M %>% lapply( function(i)
    MH_High_Dim_para_future(niter, phi[[i]], sd = c(.028, .036, .021), loglik.phi, eta[[i]], Phih, cores = 1 ))
  
  list(eta = eta, phi = phi)
}

maxi <- function(S)
{
  list(sigma2 = S%a%1,
       rho2 =   S%a%2,
       mu =     S%a%4,
       omega2 = S%a%3 - (S%a%4)^2 )
}
sigma2 rho2 mu1 mu2 mu3 omega21 omega22 omega23
Oracle 0.0499498 0.1067949 4.862881 4.035071 0.5257318 0.4788974 0.0965788 0.0126203
Initialisation 0.4219502 0.6479511 5.454432 4.454432 0.9544318 0.6285490 0.2285490 0.1385490 
niter <- 300
MH.iter <- 10
res <- SAEM(niter, MH.iter, para, Phi, S, Z, sim, maxi, burnin = 150, coef.burnin = 3/4, eps = 1e-3, verbatim = 2)
saveRDS(res, 'saem.rds')
gg <- plot(res)
## [1] "SAEM execution time = 05min 27sec"
Résultat de l’algo SAEM-MCMC
sigma2 rho2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3
Valeur réelle 0.0500 0.1000 5.0000 4.0000 0.5000 0.5000 0.1000 0.0100
Valeur estimée 0.0497 0.1358 4.9092 4.0610 0.5349 0.4602 0.0948 0.0127
Rrmse 0.0058 0.3584 0.0182 0.0153 0.0698 0.0796 0.0515 0.2655